#include "quantum.h"
#include "alloy.h"

/* --------- PHYS_VAR ---------------*/

class Phys_var {

 protected:

  double phi;  // Potential  (Volt)
  double phi_n;  // Quasi-Fermi potential for Electrons (Volt)
  double phi_p;  // Quasi-Fermi potential for Holes (Volt)
  double qphin;  // Quasi-Fermi potential for Electrons in the quantum well (Volt)
  double qphip; // Quasi-Fermi potential for Holes in the quantum well (Volt)
  double D;  // Displacemet field (C/cm^2)
  double Jn;  // Density of electron current (A/cm^3)
  double Jp;  // Density of hole current (A/cm^3)
  double alloy_x; // Fraction of the first alloy material
  double alloy_y; // Fraction of the second alloy material
  double dx;  // growth direction mesh step size
  double dy;  // in-plane direction mesh step size 
  double dt;  // time mesh step size
  double doping; // Doping density (cm^-3)
  double a;  // In plane lattice constant
  double c; // Growth direction lattice constant
	
 public:

  Alloy_var * elem; // Pointer to class containing all the material parameters
  Qntm_var * qptr_e; // Pointer to the electron quantum structure
  Qntm_var * qptr_h; // Pointer to the hole quantum structure

  // Constructors and operators
  Phys_var();
  Phys_var( const Phys_var &);  // Copy constructor
  Phys_var operator=(const Phys_var  &);

  // Dimension-related functions
  double xsize() { return dx; }
  void incr_xsize(double d) {dx = dx + d;}
  void double_xsize() {dx = 2*dx;}
  void half_xsize() {dx = 0.5*dx;}
  double ysize() { return dy; }
  double tsize() { return dt; }

  // Protected member returning functions
  double ret_alloy_x() { return alloy_x; }
  double ret_alloy_y() { return alloy_y; }
  double ret_phi() { return phi; }
  double ret_D() { return D; }
  double ret_phi_n() { return phi_n; }
  double ret_phi_p() { return phi_p; }
  double ret_qphin() { return qphin; }
  double ret_qphip() { return qphip; }
  double ret_Jn() { return Jn; }
  double ret_Jp() { return Jp; }
  double ret_a() { return a; }
  double ret_c() {return c; }
  double ret_doping() { return doping; }

  // Initialisation function
  void set ( double dx,double dy,double dt,
	     double alloy_x, double alloy_y,
	     double doping,double bias, 
	     Alloy_var * elem, Qntm_var * qptr_e, Qntm_var * qptr_h );

  // These functions modify the value of protected members
  void set_potential ( double );
  void set_a(double a_inp) { a=a_inp; }
  void set_c(double c_inp) { c=c_inp; }
  void set_c(double a_inp, double x_inp, double y_inp);
  void set_phi(double phi_inp) { phi=phi_inp; }
  void set_D(double D_inp) { D=D_inp; }
  void set_phi_n(double phi_n_inp) { phi_n=phi_n_inp; }
  void set_phi_p(double phi_p_inp) { phi_p=phi_p_inp; }
  void set_qphin(double qphin_inp) { qphin=qphin_inp; }
  void set_qphip(double qphip_inp) { qphip=qphip_inp; }
  void set_Jn(double Jn_inp) { Jn = Jn_inp; }
  void set_Jp(double Jp_inp) { Jp = Jp_inp; }

  // Functions that calculates important physical quantities
  double electric_field() {return (D-Polarization())/(EPS0*elem->eps(alloy_x,alloy_y));}
  double Polarization() {return elem->Pz(alloy_x,alloy_y,a,c)+elem->Psp(alloy_x,alloy_y);}
  double nqw();
  double pqw();
  double eps() { return elem->eps(alloy_x,alloy_y); }
  double me() { return elem->me(alloy_x,alloy_y); }
  double mh() { return elem->mh(alloy_x,alloy_y); }
  double chi() { return elem->chi(alloy_x,alloy_y); }
  double Eg() { return elem->Eg(alloy_x,alloy_y); }
  double Nc() { return elem->Nc(alloy_x,alloy_y); }
  double Nv() { return elem->Nv(alloy_x,alloy_y); }
  double En() { return -phi_n; }
  double Ep() { return -phi_p; }
  double mob_n() { return elem->mob_n(alloy_x,alloy_y); }
  double mob_p() { return elem->mob_p(alloy_x,alloy_y); }
  double Dn() { return elem->Dn(alloy_x,alloy_y); }
  double Dp() { return elem->Dp(alloy_x,alloy_y); }
  double Tn() { return elem->Tn(alloy_x,alloy_y); }
  double Tp() { return elem->Tp(alloy_x,alloy_y); }
  double cb() { return -phi-chi(); } 
  double cb(double x,double y) { return -phi-elem->chi(x,y); } 
  double vb() { return -phi-chi()-Eg(); }
  double ni() { return sqrt(Nc()*Nv())*exp(-0.5*Eg()/VT); }
  double N3D(double E) { // Return the density of states at energy E above the conduction band
    return (E>0)?(1e-6*Q)*(sqrt(Q)/PLANKUT)*me()*(MASS/PLANKUT)
      *(sqrt(MASS)/(PI*PI*PLANKUT))*sqrt(2.0*me()*E):0;
  }
  double n() { return Nc()*FD((En()-cb())/VT);}
  double n(double E, double Ec, double Ef) { return N3D(E-Ec)*GX((E-Ef)/VT,1.0);}
  double n_sum(double Ein, double Efin, double Ec, double Ef);
  double n1() { return Nc()*FD(-elem->Etrap_e()/VT);}
  double pn1() { return Nv()*FD((elem->Etrap_e()-Eg())/VT);}
  double np1() { return Nc()*FD((elem->Etrap_h()-Eg())/VT);}
  double p1() { return Nv()*FD(-elem->Etrap_h()/VT);}
  double ntrap() { return elem->Ntrap()*GX((cb()-elem->Etrap_e()-En())/VT,1.0); }
  double ntrap(double pot) { return elem->Ntrap()*GX((-pot-chi()-elem->Etrap_e()-En())/VT,1.0 );}
  double dntrap(double pot) { return -elem->Ntrap()*dGX((-pot-chi()-elem->Etrap_e()-En())/VT,1.0)/VT; }
  double ntrapf(double pn) { return elem->Ntrap()*GX((cb()-elem->Etrap_e()+pn)/VT,1.0 );}
  double dntrapf(double pn) { return elem->Ntrap()*dGX((cb()-elem->Etrap_e()+pn)/VT,1.0)/VT; }
  double ptrap() { return elem->Ptrap()*GX((Ep()-vb()-elem->Etrap_h())/VT,1.0); }
  double ptrap(double pot) { return elem->Ptrap()*GX((Ep()+pot+chi()+Eg()-elem->Etrap_h())/VT,1.0 );}
  double dptrap(double pot) { return elem->Ptrap()*dGX((Ep()+pot+chi()+Eg()-elem->Etrap_h())/VT,1.0)/VT; }
  double ptrapf(double pn) { return elem->Ptrap()*GX((-vb()-elem->Etrap_h()-pn)/VT,1.0 );}
  double dptrapf(double pn) { return -elem->Ptrap()*dGX((-vb()-elem->Etrap_h()-pn)/VT,1.0)/VT; }
  double p() { return Nv()*FD((vb()-Ep())/VT);}
  double n(double pot) { return Nc()*FD((pot+En()+chi())/VT); }
  double p(double pot) { return Nv()*FD((-pot-Ep()-chi()-Eg())/VT);}
  double nf(double pn) { return Nc()*FD((-cb()-pn)/VT);}
  double pf(double pp) { return Nv()*FD((vb()+pp)/VT); }
  double dn(double pot) { return Nc()*dFD((pot-phi_n+chi())/VT)/VT; }
  double dp(double pot) { return -Nv()*dFD((-pot+phi_p-chi()-Eg())/VT)/VT;}
  double dnf(double pn) { return -Nc()*dFD((-cb()-pn)/VT)/VT; }
  double dpf(double pp) { return Nv()*dFD((vb()+pp)/VT)/VT;}

  // Doping-related
  double ND() { return (doping>0)? doping*GD((-cb()+En()+elem->Eion(doping))/VT) : 0.0;}
  double NA() {return (doping>0)? 0.0 :-doping*GA((-vb()+Ep()-elem->Eion(doping))/VT);}
  double ND(double pot) {return (doping>0)? doping*GD((pot-phi_n+chi()+elem->Eion(doping))/VT) : 0.0;}
  double NA(double pot) {return (doping>0)? 0.0 :-doping*GA((pot-phi_p+chi()+Eg()-elem->Eion(doping))/VT);}
  double NDf(double fermi) {return  (doping>0)? doping*GD((-cb()-fermi+elem->Eion(doping))/VT) :0.0;}
  double NAf(double fermi) {return (doping>0)? 0.0 :-doping*GA((-vb()-fermi-elem->Eion(doping))/VT);}
  double dND(double pot) {return (doping>0)? doping*dGD((pot-phi_n+chi()+elem->Eion(doping))/VT)/VT :0.0;}
  double dNA(double pot) {return (doping>0)? 0.0 :-doping*dGA((pot-phi_p+chi()+Eg()-elem->Eion(doping))/VT)/VT;}
  double dNDf(double fermi) {return (doping>0)? -doping*dGD((-cb()-fermi+elem->Eion(doping))/VT)/VT :0.0;}
  double dNAf(double fermi) {return (doping>0)? 0.0 :doping*dGA((-vb()-fermi-elem->Eion(doping))/VT)/VT;}

  // Recombination functions
  double SRH_Rec() { 
    return (n()*p()-ni()*ni())/(elem->Tn(alloy_x,alloy_y)*(p()+ni())+elem->Tp(alloy_x,alloy_y)*(n()+ni()));
  }
  double SRH_Rec_ef(double imref) { 
    return (nf(imref)*p()-ni()*ni())/(elem->Tn(alloy_x,alloy_y)*(p()+ni())+elem->Tp(alloy_x,alloy_y)*(nf(imref)+ni())); 
  }
  double dSRH_Rec_ef(double imref) { 
    return (p()+elem->Tp(alloy_x,alloy_y)*SRH_Rec_ef(imref))
	    *(1.0/(elem->Tn(alloy_x,alloy_y)*(p()+ni())+elem->Tp(alloy_x,alloy_y)*(nf(imref)+ni())))*dnf(imref); 
  }
  double SRH_Rec_hf(double imref) {
    return (pf(imref)*n()-ni()*ni())/(elem->Tn(alloy_x,alloy_y)*(pf(imref)+ni())+elem->Tp(alloy_x,alloy_y)*(n()+ni())); 
  }
  double dSRH_Rec_hf(double imref) { 
    return (n()+elem->Tn(alloy_x,alloy_y)*SRH_Rec_hf(imref))
	   *(1.0/(elem->Tn(alloy_x,alloy_y)*(pf(imref)+ni())+elem->Tp(alloy_x,alloy_y)*(n()+ni())))*dpf(imref); 
  }
  double Trap_Rec_e() {
    if(elem->Tntrap()!=0 || elem->Tptrap()!=0 ) return (n()*p()-ni()*ni())/(elem->Tntrap()*(p()+pn1())+elem->Tptrap()*(n()+n1()));
    else return 0.0;
  }
  double Trap_Rec_h() { 
    if(elem->Tntrap()!=0 || elem->Tptrap()!=0 ) return (n()*p()-ni()*ni())/(elem->Tntrap()*(p()+p1())+elem->Tptrap()*(n()+np1()));
    else return 0.0;
  }
  double Opt_Rec();
  double QW_Rec();

  // The following functions are used to dump in a file all the interesting physical variables of the mesh point.
  void show_band(FILE * fp) {
    fprintf(fp,"%le  %le  %le  %le  %le  %le",cb(),vb(),-phi_n,-phi_p,cb()-elem->Etrap_e(),vb()+elem->Etrap_h());
  }
  void show_carrier(FILE * fp) {
    fprintf(fp,"%le  %le  %le  %le  %le %le  %le"
	    ,n(),p(),ntrap(),ptrap(),nqw(),pqw(),n_sum(cb(0,0),cb(0,0)+2.0,cb(0,0),-phi_n));
  }
  void show_current(FILE * fp) {
    fprintf(fp,"%le  %le  %le  %le  %le  %le  %le"
	    ,Jn,Jp,Jp+Jn,SRH_Rec(),Trap_Rec_e(),Trap_Rec_h(),Opt_Rec());
  }
  void show_potential(FILE * fp) { 
    fprintf(fp,"%le  %le  %le  %le  %le  %le  %le  %le"
	    ,phi,D,electric_field(),a,c,Polarization(),elem->Pz(alloy_x,alloy_y,a,c),elem->Psp(alloy_x,alloy_y));
  }
  void show_lev_e(FILE * fp);
  void show_lev_h(FILE * fp);
  void show_well_e(FILE * fp);
  void show_well_h(FILE * fp);
  void show_optical_bandgap(FILE * fp);


};
